################################################################################
########################17. Models using only W. Europe#########################
################################################################################

pols <- readRDS('pols_bjps')

#####Fixed Effect Model#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)

west1 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc1 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc1 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc1: scale_socioeconomicsal +
               scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols[pols$eastwest == 'west',], 
             na.action = na.omit, 
             model = "within")

coeftest(west1, vcov = function(x) vcovBK(x, cluster="group"))

west2 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc2: scale_socioeconomicsal +
               scale_ggini_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols[pols$eastwest == 'west',], 
             na.action = na.omit, 
             model = "within")

coeftest(west2, vcov = function(x) vcovBK(x, cluster="group"))

west3 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc3: scale_socioeconomicsal +
               scale_ggini_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols[pols$eastwest == 'west',], 
             na.action = na.omit, 
             model = "within")

coeftest(west3, vcov = function(x) vcovBK(x, cluster="group"))

west4 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc4: scale_socioeconomicsal +
               scale_ggini_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols[pols$eastwest == 'west',], 
             na.action = na.omit, 
             model = "within")

coeftest(west4, vcov = function(x) vcovBK(x, cluster="group"))

west5 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc5: scale_socioeconomicsal +
               scale_ggini_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols[pols$eastwest == 'west',], 
             na.action = na.omit, 
             model = "within")

coeftest(west5, vcov = function(x) vcovBK(x, cluster="group"))

##### Only West European Cases: Point estimates#####

fu1 <- as.data.frame(summary(west1, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(west2, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(west3, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(west4, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(west5, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 83)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#Model Fit Statistics#

fit1f <- c(summary(west1)$r.squared, summary(west1)$fstatistic$statistic, summary(west1)$fstatistic$p.value)
fit2f <- c(summary(west2)$r.squared, summary(west2)$fstatistic$statistic, summary(west2)$fstatistic$p.value)
fit3f <- c(summary(west3)$r.squared, summary(west3)$fstatistic$statistic, summary(west3)$fstatistic$p.value)
fit4f <- c(summary(west4)$r.squared, summary(west4)$fstatistic$statistic, summary(west4)$fstatistic$p.value)
fit5f <- c(summary(west5)$r.squared, summary(west5)$fstatistic$statistic, summary(west5)$fstatistic$p.value)


fit.statsf <- rbind(fit1f, fit2f, fit3f, fit4f, fit5f)

print(mean(fit.statsf[,1])) #R-Squared
print(mean(fit.statsf[,2])) #Adjusted R-Squred
print(mean(fit.statsf[,3])) #F-Statistic
print(mean(fit.statsf[,4])) #p-value of F-Statistic

test <- as.data.frame(mvrnorm(1000, mu =coef(west1), Sigma =vcovHC(west1)))

pols_rug <- readRDS('pols_bjps') %>%
  filter(eastwest == 'west') %>%
  dplyr::select(scale_econsdw,scale_dispgini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_ggini_inc1[pols$eastwest == 'west'], na.rm = T), to = max(pols$scale_ggini_inc1[pols$eastwest == 'west'], na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*-1 + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.west <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = -1", x = "", y = "") +
  theme_bw ()+
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*0 + test[j,10]*marg[i,1]*0
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.west <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*1 + test[j,10]*marg[i,1]*1
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.west <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

fe_west_euro_sortx <- grid.arrange(lwr.plot.west, middle.plot.west, righ.plot.west, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Sorting by Income")

#ggsave("figure_A7.pdf", plot = fe_west_euro_sortx, width = 11, height = 5, units = 'in')

######Figure 5: Salience on X with only west european cases#####################

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_socioeconomicsal[pols$eastwest == 'west'], na.rm = T), to = max(pols$scale_socioeconomicsal[pols$eastwest == 'west'], na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*-1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sal <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.w <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*0 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*0
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.w <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*1
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(al = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.w <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 4.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x =scale_socioeconomicsal), inherit.aes = F)

fe_west_euro_salx <- grid.arrange(lwr.plot.w, middle.plot.w, righ.plot.w, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Salience of Economic Issues")

#ggsave("figure_A8.pdf", plot = fe_west_euro_salx, width = 11, height = 5, units = 'in')





#####Multi-Level Models - Not in text or supplementary materials#####

library(nlme)

#####Only Western European Countries#####

west1 <- lme(scale_econsdw ~ scale_ggini_inc1 +
               scale_socioeconomicsal +
               scale_ggini_inc1 *scale_socioeconomicsal +
               scale_dispgini +
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp +
               scale_ggini_inc1 *scale_dispgini +
               scale_socioeconomicsal * scale_dispgini +
               scale_ggini_inc1 * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 |  country.name , 
             data = pols[pols$eastwest == 'west',],
             na.action = na.omit,
             method = 'ML')

west2 <- lme(scale_econsdw ~ scale_ggini_inc2 +
               scale_socioeconomicsal +
               scale_ggini_inc2 * scale_socioeconomicsal +
               scale_dispgini +
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp +
               scale_ggini_inc2 * scale_dispgini +
               scale_socioeconomicsal * scale_dispgini +
               scale_ggini_inc2 * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 | country.name , 
             data = pols[pols$eastwest == 'west',],
             na.action = na.omit,
             method = 'ML')

west3 <- lme(scale_econsdw ~ scale_ggini_inc3 +
               scale_socioeconomicsal +
               scale_ggini_inc3 * scale_socioeconomicsal +
               scale_dispgini +
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp +
               scale_ggini_inc3 * scale_dispgini +
               scale_socioeconomicsal * scale_dispgini +
               scale_ggini_inc3 * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 |  country.name , 
             data = pols[pols$eastwest == 'west',],
             na.action = na.omit,
             method = 'ML')

west4 <- lme(scale_econsdw ~ scale_ggini_inc4 +
               scale_socioeconomicsal +
               scale_ggini_inc4 * scale_socioeconomicsal +
               scale_dispgini +
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp+
               scale_ggini_inc4 * scale_dispgini +
               scale_socioeconomicsal * scale_dispgini +
               scale_ggini_inc4 * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 |  country.name , 
             data = pols[pols$eastwest == 'west',],
             na.action = na.omit,
             method = 'ML')

west5 <- lme(scale_econsdw ~ scale_ggini_inc5 +
               scale_socioeconomicsal +
               scale_ggini_inc5 * scale_socioeconomicsal +
               scale_dispgini +
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp +
               scale_ggini_inc5 * scale_dispgini +
               scale_socioeconomicsal * scale_dispgini +
               scale_ggini_inc5 * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 | country.name , 
             data = pols[pols$eastwest == 'west',],
             na.action = na.omit,
             method = 'ML')

#####Point Estimates Using only West European Countries#####

r1o <- r.squaredGLMM(west1)
r2o <- r.squaredGLMM(west2)
r3o <- r.squaredGLMM(west3)
r4o <- r.squaredGLMM(west4)
r5o <- r.squaredGLMM(west5)

r.o <- rbind(r1o, r2o, r3o, r4o, r5o)
print(mean(r.o[,1]))
print(mean(r.o[,2]))

fit1 <- c(summary(west1)$logLik, summary(west1)$AIC, summary(west1)$BIC)
fit2 <- c(summary(west2)$logLik, summary(west2)$AIC, summary(west2)$BIC)
fit3 <- c(summary(west3)$logLik, summary(west3)$AIC, summary(west3)$BIC)
fit4 <- c(summary(west4)$logLik, summary(west4)$AIC, summary(west4)$BIC)
fit5 <- c(summary(west5)$logLik, summary(west5)$AIC, summary(west5)$BIC)

fit.stats <- rbind(fit1, fit2, fit3, fit4, fit5)

print(mean(fit.stats[,1])) #log liklihood
print(mean(fit.stats[,2])) #AIC
print(mean(fit.stats[,3])) #BIC

#Standard deviation of random intercept from summaries above#

print(mean(c(0.294041, 0.3005514,  0.3003577, 0.2997737, 0.2904789)))#country

#####Point Estimates#####
m1o <- as.data.frame(summary(west1)$tTable) 

m2o <- as.data.frame(summary(west2)$tTable) %>%
  dplyr::rename(Value2 = Value, 
                Std.Error2 = Std.Error, 
                DF2 = DF, 
                `t-value2` = `t-value`, 
                `p-value2` = `p-value`)
m3o <- as.data.frame(summary(west3)$tTable) %>%
  dplyr::rename(Value3 = Value, 
                Std.Error3 = Std.Error, 
                DF3 = DF, 
                `t-value3` = `t-value`, 
                `p-value3` = `p-value`)
m4o <- as.data.frame(summary(west4)$tTable) %>%
  dplyr::rename(Value4 = Value, 
                Std.Error4 = Std.Error, 
                DF4 = DF, 
                `t-value4` = `t-value`, 
                `p-value4` = `p-value`)
m5o <- as.data.frame(summary(west5)$tTable) %>%
  dplyr::rename(Value5 = Value, 
                Std.Error5 = Std.Error, 
                DF5 = DF, 
                `t-value5` = `t-value`, 
                `p-value5` = `p-value`)

dot.data.o <- cbind(m1o,m2o, m3o,m4o,m5o)

dot.data.o <- dot.data.o[ , order(names(dot.data.o))]

dot.table.o <- dot.data.o %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.o, Value, Value2, Value3, Value4, Value5)),
            se = rowMeans(dplyr::select(dot.data.o, Std.Error, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            DF = rowMeans(dplyr::select(dot.data.o, DF, DF2, DF3, DF4, DF5)),
            t.value = rowMeans(dplyr::select(dot.data.o, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.o, `p-value`, `p-value2`, `p-value3`, `p-value4`, `p-value5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.o$Predictor <- c("Intercept", "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)", "ENEP", "GALTAN Polarization", "District Magnitude", "Unemployment", "Income.Dif:Salience", "Income.Dif:Gini", "Salience:Gini", "Income.Dif:Salience:Gini")

dot.plot.o <- dot.table.o

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, 
                               levels = c( "Income.Dif", "Socio-Economic Salience", "Disposable Income Inequality (Gini)",  "GALTAN Polarization","ENEP", "District Magnitude", "Unemployment", "Income.Dif:Salience", "Income.Dif:Gini", "Salience:Gini", "Income.Dif:Salience:Gini", "Intercept"))

dot.plot.o$Predictor <- factor(dot.plot.o$Predictor, levels=rev(levels(dot.plot.o$Predictor)))

dot.plot.o %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 78)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table.o) <- dot.table.o$Predictor
dot.table.o <- dot.table.o[,-8]

print(xtable::xtable(dot.table.o, digits = 3))